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1  J  I 


We  present  an  algorithm  for  the  evaluation  of  the  Fourier  transform  of  piecewise  constant 
functions  of  two  variables.  The  algorithm  overcomes  the  accuracy  problems  associated  with 
computing  the  Fourier  transform  of  discontinuous  functions;  in  fact,  its  time  complexity 
is  0(log^{l/c)A'^  +  log^(l/c)jV^log  jV),  where  c  is  the  accuracy  and  N  is  the  size  of  the 
problem.  The  algorithm  is  based  on  the  Lagrange  interpolation  formula  and  the  Green's 
theorem,  which  are  used  to  preprocess  the  data  before  applying  the  Fast  Fourier  transform. 
It  admits  natural  generalizations  to  higher  dimensions  and  to  piecewise  smooth  functions. 
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1 


Introduction 


The  Fast  Fourier  Transform  (see,  for  example,  [2,  3,  8])  is  a  ubiquitous  tool  of  numerical 
analysis,  essential  in  signal  processing,  electrical  engineering,  VLSI  circuit  modeling,  medical 
imaging,  and  innumerable  other  fields.  It  is  a  very  effective  tool  when  the  function  to  be 
transformed  is  smooth:  however,  if  the  function  has  a  jump  discontinuity,  then  the  accuracy 
of  the  FFT  is  significantly  reduced.  This  loss  of  accuracy  is  the  result  of  the  simple  fact 
that  FFT,  from  the  mathematical  point  of  view,  is  a  collection  of  integrals  evaluated  with 
trapezoidal  rule.  Therefore,  the  numerical  error  that  results  from  the  discontinuity  is  on  the 
order  of  n~^  for  a  problem  of  size  n.  Such  errors  make  even  single  precision  calculations, 
especially  in  higher  dimensions,  prohibitively  costly.  The  cost  can  be  reduced  with  the 
help  of  the  Richardson  extrapolation,  but  double  precision  calculations  are  still  virtually 
impossible. 

We  construct  an  algorithm  that  successfully  attacks  this  problem;  for  example,  on  a 
■512  X  51?  array  our  implementation  takes  about  .50  times  the  CPU  time  of  one  512  x  512 
FFT  for  single  precision  results  and  about  160  times  for  double  precision. 

The  algorithm  uses  Green's  theorem  to  replace  the  area  integrals,  that  naturally  occur  in 
computing  the  two-dimensional  Fourier  transform  of  piecewise  constant  functions,  with  line 
integrals,  thereby  reducing  significantly  the  number  of  nodes  required  for  integration;  after 
that  the  weights  are  redistributed  onto  the  uniform  grid  with  the  help  of  the  Lagrange  inter¬ 
polation  formula,  bringing  the  data  to  the  form  suitable  for  the  FFT.  A  detailed  description 
of  the  algorithm  can  be  found  in  Section  4.  and  performance  results  in  Section  5. 

In  this  paper,  we  present  the  analysis  in  two  dimensions  under  the  assumption  that  the 
functions  to  be  transformed  are  piecewise  constant  and  that  the  discontinuities  occur  along 
a  collection  of  polygons.  Neither  of  these  is  a  serious  limitation;  the  algorithm  generalizes 
quite  naturally  to  higher  dimensions  and  to  curved  boundaries.  In  Section  6  we  describe  a 
generalization  of  our  algorithm  to  piecewise  smooth  functions. 

Acknowledgements:  The  author  would  like  to  thank  Professors  Ronald  Coifman  and 
Vladimir  Rokhlin  for  several  useful  discussions.  He  would  al-so  like  to  thank  Professor  Eytan 
Barouch  for  suggesting  the  problem  and  supplying  the  VLSI  masks  for  the  experiments. 

2  Mathematical  Preliminaries 

2.1  Lagrange  Interpolation 

The  Lagrange  interpolation  formula  (see,  for  example,  [1,  4])  approximates  a  function  /  : 
]R’  —  f '  by  the  expression; 

P  P 

fix)  =  £  fiXrr,)  n  (') 

m  =  l 

n#m 


1 


where  ji  <  . . .  <  Xp  are  points  on  the  real  line  and  Rpix)  is  the  error  term.  For  each 
m  =  1,. .  .,p,  we  wiU  denote  by  6^  the  polynomial  defined  by  the  formula: 


and  observe  that 


6m(x)=  n 


X  -  x„ 


m 


^m{Xn)  -  I  0  i 


if  m  =  n, 
if  m  ^  n. 


In  view  of  (2),  formula  (1)  can  be  rewritten  as 


(2) 


f{x)=  ^  +  Rp{x) 


Tn=l 


(3) 


Rp[x)  in  equations  (1)  and  (3)  is  the  error  term  and 

f^p(x)  =  - j —  for  some  5€(xi,Xp). 

P.  ^  * 

n=l 

In  our  case,  the  nodes  xj,...,Xp  are  going  to  be  uniformly  spared  with  the  sampling 
intervd  /i  >  0.  and  the  point  x  will  lie  in  the  interval  of  length  h  centered  on  the  center  of 
the  interpolation  window  [xi,Xp].  Furthermore,  the  function  /  will  be  of  the  form  /(x)  = 
exp(2jrifcx).  Under  these  conditions. 


l«p(x)l  <  { 


1  .  lI-iil-LlL  .  for  odd  p. 

W  Tl' 


[r(f  Fi))‘ 

p! 

(r(f +  1)]' 


h^{2itk)^  for  even  p 


h^{2Trkf, 


where  r(r )  is  the  Gamma  function:  Ff ;  +  1 )  =  r!. 

We  will  also  make  use  of  the  two-dimensional  Lagrange  interpolation  formula: 


p  p 


(>m(x)  ■  ■  f(Xn,yr,)  +  Rp{x,y). 


n=l  m  =  l 


(•'>) 


Here  the  points  {(xm-yn)}m.n=i  interpolation  nodes,  the  functions  /ij . /ip  are  the 

same  as  in  (2),  and  the  error  term.  Rp{x.y),  admits  the  bound 


\Rp(x,y)\<  ■  [h;(2>^<:)»-^hg(27^1^] 


(6) 


when  the  function  /(x,y)has  the  form  exp(27r*(A:x-l-/y))  and  the  nodes  are  uniformly  spared 
with  the  sampling  interval  hj  along  the  x-axis  and  hy  along  the  y-axis. 
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2.2  Green’s  theorem 

We  will  encounter  area  integrals  of  the  form 


L 


^-2irimx  ^-Irtny 


dydi. 


where  D  is  a  bounded  domain  in  R^.  We  will  replace  them  with  line  integrals  using  the 
following  version  of  Green's  theorem  in  the  plane: 

Here  Q  :  R^  — >  €*  is  a  function  in  C^(D)  and  F  is  the  boundary  of  D  traversed  counter 
clockwise.  Applying  (7)  to  Q(x,y)  —  we  arrive  at  a  formula  that  will  be 

important  to  us  later: 


Jd 


_J_  f  ^-2’r,mx^-2mny^  ^  ^  q 

-27rim  Jr  »  r 

^  a-c-2’'*"V  dy 


when  m  =  0. 


2.3  Gauss-Legendre  Quadratures 

Gauss- Legendre  integration  scheme  (see.  for  example.  [4,  7])  provides  an  approximation  to 
the  integral  of  a  function  /  :  [0, 1]  —  O’: 

/  /(t)df=  /(t*)^k+ £,{/).  (9) 

*=1 


where  the  points  t\ . t,  are  the  Gauss-Legendre  nodes,  the  numbers  . .j,  are  the 

Gauss-Legendre  weights,  and  £,(/)  is  the  error  term.  The  Gauss-Legendre  nodes  and 
weights  are  chosen  to  make  the  approximation  (9)  exact  whenever  the  function  /  is  a 
polynomial  of  degree  less  than  2q.  Consequently,  the  Gauss-Legendre  scheme  is  effective 
for  functions  that  are  well  approximated  by  polynomials;  in  fact: 


£,(/)  = 


(2<?-L  l)[{2v)!P 


The  proof  of  ( 10)  can  be  found  in  [4]. 


for  some  ^  £  [0.  1]. 


(10) 


3  Mathematical  Description  of  the  Algorithm 

3.1  Statement  of  the  Problem 

Our  goal  in  this  paper  is  to  find  a  way  to  compute  Fourier  Transforms  of  piecewise  constant 
functions  in  R^.  In  other  words,  we  would  like  to  construct  an  algorithm  that  for  every 


.3 


piecewise  constant  function  /  :  — *  C'  will  compute  an  approximation  to  the  collection 

of  integrals,  to  which  we  wiD  refer  as  the  Fourier  Transform  /  of  /: 

/(m,n)1l^  r  (11) 

Jo  Jo 

where  m  and  n  are  integers  such  that  -M  <  m  <  M  and  -N  <  n  <  N .  We  will  refer  to 
the  pair  (m.n)  as  the  frequency.  The  pair  (M,N),  thus,  is  the  highest  frequency  of  interest. 
The  most  common  approximation  to  the  integrals  (11),  the  Discrete  Fourier  Transform  [2]  is, 
essentially,  the  trapezoidal  rule;  its  standard  implementation  is  the  Fast  Fourier  Transform 
(see,  for  example,  [2,  3,  8]).  For  discontinuous  functions  the  trapezoidal  rule  has  accuracy 
on  the  order  of  for  an  ,V  by  N  problem,  or,  conversely,  for  accuracy  e,  will  require 
on  the  order  of  operations;  such  performance  is  woefully  inadequate  in  most  practical 
problems. 

In  this  paper  we  solve  the  following  problem: 

Problem  1  For  every  piecewise  constant  function  f  with  discontinuities  along  a  finite  num¬ 
ber  of  polygons,  and  for  every  accuracy  e  >  0,  compute  all  the  integrals  (I  I)  with  accuracy 
f  in  the  number  of  operation  not  greater  than 


O  ^  ^log^  .V  ^  +  ^log^  7)  ^  'V  ^  • 


where  S  =  max{.V/,.V}. 

.Mathematical  description  and  error  analysis  of  the  algorithm  that  solves  Problem  1  can  be 
found  in  §3.2  and  §3.3,  respectively.  Section  4  contains  the  full  description  of  the  algorithm, 
and  Section  5  contains  the  results  of  some  of  our  numerical  experiments.  In  Section  6  we 
describe  an  extension  of  the  algorithm  to  piecewise  smooth  functions. 

3.2  Description  of  the  Algorithm 

A  piecewise-constant  function  /  :  E^  ->  €'  can  be  written  as  a  Unear  combination  of 


characteristic  functions.  By  the  characteristic  function  I5  :  E^ 
mean,  as  usual: 


C’  of  a  set  .S  C  E^  we 


Using  this  notation,  we  write  the  piecewise-constant  /  as: 


^  A'j  •  l/5^(x.y).  (13) 

r=l 

where  J  is  a  positive  integer  constant  determined  by  the  problem,  and,  for  each  j,  Kj  is 
a  complex  constant  and  Dj  is  a  bounded  domain  in  E^.  We  wiU  assume  from  now  on 
that  all  domains  Dj  Ue  within  the  unit  square,  their  interiors  are  pairwise  disjoint,  and  the 
boundaries  r_,  of  the  domains  Dj  are  polygons.  Indeed,  for  each  j,  let  {7/  :  [0. 1]  — 
be  the  set  of  paths  that  form  F^.  Then,  our  assumptions  are: 


4 


t 


1-  Oj  C  [0, 1]  X  [0, 1]  for  all  j, 

2.  Int(Dj)  n  lnt{Dk)  =  0  for  till  j  ^  k. 

3.  Every  path  7/  has  the  form: 

=  {ao  +  at,bo  +  bt)  for  t  €  [0, 1].  (14) 


The  first  two  assumptions  do  not  restrict  generality  at  all,  and  the  third  one  is  made  to 
simplify  the  error  estimates  and  the  implementation.  It  does  not  affect  the  algorithm  itself. 
Substituting  (13)  into  (11)  and  using  (8),  we  get: 

f{m.  n)  =  /’  /'  E  ■  1d,(x,  dydx 

Jo  Jo  ^ 

J  f 

=  E  /  dydx 

=  E/  K:  ■  F'm.n{x,y)dy,  (15) 

j=i  ■'^1 


where  each  is  the  boundary  of  the  domain  Dj  and  the  function  f „  „  :  IR^ 
defined  by  the  formula: 


/  ^-2)rimr^-2)riny 


F  i 

•  m,n  —  \ 


X€ 


— 27rim 

-2irtnv 


when  m  7:  0. 
when  m  =  0. 


r'  is 

(16) 


For  each  j.  let  {7/}  be  the  set  of  paths  that  form  Fj.  Along  each  7/  we  approximate  the 
line  integrals  in  (15)  with  the  help  of  Gauss-Legendre  quadrature  scheme  of  '  2.3.  Indeed, 
using  approximation  (9)  on  each  7/  and  substituting  the  result  into  (15)  we  get: 

f(m,n)  =  EEE  ■  b  ■  ujI.  (17) 

j=i  (=1  fc=i 

where,  for  each  /,  the  points  {( V’*)}*-,  are  the  images  under  7/  of  the  Gauss-Legendre 
nodes  tk  of  equation  (9): 

=  llUkh  (18) 

the  numbers  are  the  corresponding  Gauss-Legendre  weights,  and  6  has  been 

defined  in  equation  (14). 

In  order  to  be  able  to  use  the  Fast  Fourier  Transform,  we  redistribute  the  weights 
from  the  Gauss-Legendre  nodes  in  (18)  to  the  uniform  grid  with  the  help  of  the  Lagrange 
Interpolation  formula  (5)  applied  to  the  function  Fjn,n  at  (x,y)  =  (x*,  V’*)  for  each  k  and  /. 
The  final  approximation  to  f(m,n)  becomes: 

/(m.n)  w  5IEE  ^'2  •<’•4  E  E  (*9) 

j=i  1=1  t=i  ;i=ij2=i 


♦ 


I 


» 


» 


» 


5 


where,  for  each  k  and  /,  the  points  ^te  chosen  on  the  uniform  grid  so  that 

the  point  (\[.  t'’[)  lies  in  the  hj-  by  hy  rectangle  centered  on  the  center  of  the  interpolation 
window  [xt.x*]  X 

Remark  3.1  Changing  the  order  of  summation  in  (19),  we  see  that 

(J  71  \ 

J=1  t=l  k=l  / 

where  the  points  {(Xj, .  Jtjj )}  form  a  uniform  grid  on  [0, 1]  x  [0, 1]  with  sampling  intervals  hj 

and  hy. 

Let  G  :  H?  — *  C’  be  defined  by  the  term  in  the  parentheses  in  (20): 

G{j..  j2)  ='  E  E  E  AV  fc  •  4  ■  ^..(xl-)  •  ^22(4).  (21 ) 

j=i  /=i  t=i 

Then,  substituting  (21)  into  (20),  and,  using  the  definition  of  Tm.n  in  (16),  we  get; 

/(m.n)  «  EE  G(jl,j2)  ■  Ani,n(  Xji ,  ) 


-2irim 


- — EE G( ji. ja)  •  f  vvhen  m  ^  0. 

im 


n  \  n  ) 


when  m  =  0 


which  can  be  readily  computed  for  aU  rn  and  n  with  the  help  of  two-  and  one-dimensional 
FFTs.  ■ 

3.3  Error  Estimates 

In  this  section  we  estimate  the  accuracy  of  the  approximations  ( 17)  and  ( 19)  made  in  §3.2. 
Let  ij  be  a  straight  path  as  defined  in  (14): 

tjit)  =  {ao  +  at.bo  +  bt)  for  <  €  [0.  1]. 

Since  Gauss- Legendre  quadrature  and  Lagrange  interpolation  formula  are  translation  in¬ 
variant,  the  estimates  will  not  change  if  we  assume  that  ao  =  bo  —  0.  Then,  for  every  (m.v) 
such  that  -  A/  <  m  <  M  and  -  A’  <  n  <  .V,  using  (9)  and  (14),  we  get: 

/  F„.n(x,y)dy  =  f  FmA^^bt)  ■  bdt 
J-rl  Jo 

=  E  ^"•'‘(xlt'4) 

fc=i 

+  bEy,{Fr„,„otl).  (23) 


Putting  (10),  (28),  and  (29)  together,  we  get  an  estimate  of  the  error  of  Gauss- Legend  re 
quadrature  in  (2d): 


< 


(2?/ +  1)[(27()!P  2jr|m| 


(2qi  +  1)[(29/)!P 


■(2TMlf‘>‘  qi 


(29,  1)[(29,)!]3 


3L 

2ir' 


which  is  exactly  the  clciim  of  the  lemma.  ■ 


when  ju  ^  0, 


when  rn  ^  0 


Corollary  3.1  Let  y  :  [0,  1]  —  IR^  he  a  straight  path  defined  by  (14),  and  let  F^.n  be  a 
function  defined  by  (16). 

The  n.  for  e  very  f  >  0. 


y 

I  f^m.n  dy  —  ^  ~  l^m,n{~){lk))  '  b  ■ 


k=\ 

whenever  —M  <  ni  <  .\/,  -.V  <  r,  <  A’,  and 

q  >  const.  ■  max 


<  e 


{A/Z,logl|. 


(30) 


(31) 


Here  te  and^'k  are  the  (laus.^^- Legendre  nodes  and  weights,  respectix'ely.  and  A/,  L  hax'e  been 
defined  in  ('2*'/. 


Proof.  I'sing  Stirling's  appro.ximation  to  9!.  we  see  that  for  large  9  the  right-hand  side  of 
expression  (21)  is  asymptotic  to: 


2'‘^+'(9!)‘'7r2’9 
(29-h  l)[(29)!]^27r 


Therefore,  the  inequality  (30)  will  hold  if 


q 


>  const.  ■  max 


(32) 


(33) 


Remark  3.2  The  quantity  M  is  the  frequency  of  the  function  F\f  ,\  along  y,  and  /,  is  the 
length  of  7.  Corollary  3.1.  then,  implies  that  the  number  of  Gauss- Legendre  nodes  needed 
to  compute  the  integral  Fm  s  dy  with  accuracy  e  is  proportional  to  the  number  of  periods 
of  Lm.\  along  7  (so  long  as  that  number  is  greater  than  log(l/c)).  ■ 


.After  computing  Gauss- Legendre  nodes  and  weights,  in  the  equation  ( 19)  we  redistribute 
the  weights  onto  the  uniform  grid  with  the  help  of  Lagrange  interpolation  formula: 

p  p 

^m,n(  Xk- t'T  )  =  bjj(Xk)  ■  bj.j{xl’f.)  ■  Fm.niTj^.yj^)  +  Rp{\‘f..X.'[).  (.14) 

Ji  =  1  J2  =  I 

The  error  Rp{  V’[)  of  this  approximation  is  estimated  in  the  next  lemma. 
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•  4 


•  I 


•  4 


»  • 


»  4 


9  4 


»  4 


»  4 


»  4 


• _  •  •  •  •  • 


•  L 


Lemma  3.2  tor  alt  (rji.Ti)  such  that  -M  <  m  A/  atul  <  a  ■  .\  ,  tin  t  t  r  (tr  ft  f  jt, 

Rjt  \  .1')  satisfiis  tin  ituqualily: 

[r  (i-'  +  Ill'* 

lUpW-v)]  <  --p-,  +  ,:r. 

<  ,^/Trp  ■ 1^']  i.ihi 

for  p  >  '2. 

Proof.  For  m  ^  0.  (3o)  follows  dirertly  from  (fil.  When  in  =  0.  v.f  observe  that  th'- 
function  x  is  approximated  exactly  hy  Lai^range  interpolation  formula  for  p  2  <  Hud  we 
only  consider  those),  so  the  estimate  (H'>)  applies  to  this  case  also.  Imupialit^  i  dti  i  i-  a 
consequence  of  Stirlin&'s  fortiiula  'll.  ■ 

Definition  3.1  Let  =  I'lp.fi  he  the  minimum  numher  of  points  p<'r  w aveleiiut h  that  i- 
needed  to  make  LagraiiRe  interpolat  loti  fortiiula  of  or<ler  p  approximate  w  it  h  ar  iuraii  >  all 
the  function' 

X,  V  ,,  ^,,,1  _  y  ,  .  y 

inside  the  central  h,  hy  rect  atigle  of  t  he  interpolat  loii  window  ■  ’,vi .  V;  ■ 

Remark  3.3  In  orcler  to  assure  that  l.agratiee  interpolation  has  ac  i  urai  \  . .  we  now  i  hoo-. 
t  he  parameter' p  and  /-.defined  in  eip/ation  i  I  i  atid  Definition  d,  1 .  re'pei  nveK  In  view 
Definition  d.l.  wi-  set  the  sampline  intervals  and  h.^  to  he; 

h,  =  -TV  and  ot, 

I'M  ('.\ 

Stibstitiitine  I -IT )  into  (d(i).  we  see  tliat 

^  ■  [I'J-h^A/ I'  ♦  Ci-r/i^.V  d;  •  .^.Tp  {  -  )  i  .p> 

.Numerical  evidence  stiK^ests  that  r  =  p'2  i'  a  good  rhoiee  .Making  thi-  (fioi/e  u;  t'  we 
see  that  Lagrange  interpolation  is  guaranteed  to  have  aicuraiv  <  if 


Consequently,  we  choose  parameters  p  and  r  of  the  computation  accordinc  to  the  relatioi  - 

2i'  -  p  log  - .  I  HI 

« 

Numerical  experiments  of  Section  o.  however,  sliow  that  the  requirement  i  .ft ,  o  undulv 
strict.  For  example,  for  *  =  10"'‘*.  the  inequality  (dft)  demands  p  -  2‘>.  while  p  ;  If,  i- 
suffirient.  Similarly,  for  single  precision,  t  =  lO"'.  the  inequality  i  d'l  i  require'  p 
when,  in  fact,  p  =  10  is  large  enough.  ■ 


With  the  help  of  the  lemmas  3.1  and  3.2  we  can  now  estimate  the  total  error  of  the 
computation. 

Substituting  (34)  into  (23),  we  get; 

I  y)dy  =  H  H  .  Vj*  )  •  t  •  4 

•''^1  fc=ljl=lj2  =  l 

+  6-£„(F„,„o7/)  +  '^b  Jk  -Rpixi.y’i)  (41) 

*=1 

Denoting  by  Rj,  the  right-hand  side  of  (36),  by  £„  the  right-hand  side  of  (24 ),  and,  observing 
that  ==  1,  we  see  that  the  error  of  integrating  along  the  path  7/  in  (41)  is  bounded 

by 

6£„ -h  6ftp  <  26f  <  2L<,  (42) 

if  p  and  1/  are  chosen  according  to  relations  (40)  in  Remark  3.3. 

Summing  up  the  errors  in  (42)  over  all  paths  -yj  that  form  the  boundary  Pj.  we  see  that 
integration  around  each  Pj  results  in  an  error  not  greater  than 

2^lir,||. 


where  ||Pj|l  denotes  the  length  of  the  perimeter  of  Pj.  The  total  error  of  the  computation 
is,  therefore,  at  most 

J 

J=1 

Actual  errors  are  reported  in  Section  .5. 


4  Formal  Description  of  the  Algorithm 

4.1  Description  of  the  Algorithm 

Let  /  :  IR^  — ►  C*  be  as  in  (13),  and  M  and  N  be  the  highest  frequencies  desired  in  the 
Fourier  transform  of  /  along  the  x-  and  j/-axis,  respectively. 

Step  0:  Initialization. 

Comment:  [  We  choose  precision  of  the  computation  f  and  determine  the  degree  of 
Lagrange  interpolation  p  and  the  oversampling  factor  1/  as  described  in  Remark  3.3. 
Since  we  will  need  i>  points  per  wavelength,  we  create  the  function  G(m,n).  defined  in 
(21),  on  the  uniform  grid  on  [0, 1]  x  [0, 1]  with  the  sampling  intervals  )  and 

hy  =  lf{i/N).  We  will  also  make  use  of  a  one-dimensional  array  G’o(n)  to  compute  the 
Fourier  Transform  for  the  case  m  =  0.  We  precompute  the  denominators  of  Lagrange 
interpolation  and  Gauss- Legendre  nodes  and  weights  on  [0, 1).  ] 


t 
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do  H  =  1, . . . ,  i>N 
•  Go(n)  =  0 

do  m  =  1, . .  .,i/M 
G(m,n)  =  0 

<  enddo 

enddo 

Step  l:  Green’s  Theorem  and  Lagrange  Interpolation. 

Comment:  [  For  each  j,  we  integrate  the  constant  function  A'j  along  each  of  the  paths 
■yf  :  [0, 1]  — ►  Fj,  defined  in  (14),  and  then  redistribute  the  weight  from  each  of  the  re¬ 
sulting  Gauss- Legendre  nodes  to  Lagrange  nodes  from  the  uniform  grid  as  described 
in  equation  (19).  Note:  The  double  index  (m(ii),  n(t2))  of  the  array  G  corresponds  to 
the  point  (x„,y,J.  ] 

do  j  =  1,. .  .,J 
do  /  =  1 , . . . ,  Lj 

Using  the  precomputed  nodes  and  weights  on  [0, 1]  calculate  the  Gauss- 
Legendre  nodes  {(x!f  t^’jt)}  a^nd  weights  u;[  (defined  in  (18))  needed  to  com¬ 
pute  the  line  integrals  along  yf. 
do  k  =  1 . . . . ,  9/ 

do  Ji  =  1, _ p  and  =  1, _ p 

6'(m(ji).n(i2))  =  6’(m(ii),  n(i2))  +  A';  '  biA^'k) 

Go(ri(t2))  =  6'o(n(i2))  +  A'j  •  ■  h  ■  ■  d.jivi) 

enddo 

enddo 

enddo 

enddo 

Step  2;  Foirier  Transform. 

Comment;  [  We  apply  the  two-dimensional  FFT  to  the  i/M  by  l'.V  array  G  and  only 
keep  the  frequencies  (m.n)  with  -M  <  m  <  M  and  -A’  <  n  <  .V.  We  store  the 
result  in  the  array  FG.  Similarly,  we  apply  the  one-dimensional  FFT  to  the  array 
Go.  keep  only  the  frequencies  n  with  —.V  <  n  <  A’,  and  store  the  result  in  the  array 
FGo.  ] 

Step  3:  Integration. 

Comment:  [  For  each  frequency  (m,  n)  with  m  /  0,  we  divide  its  coefficient  by  -27rim. 
which  brings  the  computation  in  agreement  with  equation  (19).  The  case  m  =  0  is  jil- 
ready  correct  and  needs  no  adjustment.  ] 

dom=-M-fl,...,M  and  n  =  -N  -1-  1 , . . . ,  A' ;  with  m  ^  0 
FG(m,n)  =  FG{m.n)/(-2irim) 

enddo 

don  =  —N  -f  1 . N 

FG{0,n)=  FGoin) 

enddo 

II 


» 


» 
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Remark  4.1  Often,  for  example  in  VLSI  modeling,  the  boundaries  of  the  domains  Dj  con¬ 
tain  vertical  and  horizontal  straight-line  segments.  Contribution  from  horizontal  segments 
to  (15)  is  zero,  and  if  5  is  a  vertical  segment  from  (xq,  i/o)  to  (xo.yi),  then,  for  n  ^  0. 


L 


{x,y)dy 


^m,n(^0i  y\  )  ^ni,n{^0»  J/o) 

— 27rin 


(43) 


that  is.  computation  of  a  line  integral  is  replaced  by  evaluation  of  a  sum  of  two  terms  of 
the  same  type  as  in  the  approximation  (17).  Effects  of  the  resulting  speed-up  are  discussed 
in  Section  5.  ■ 


4.2  Complexity  Analysis  of  the  Algorithm 

Time  complexity  of  the  algorithm  is  summarised  in  Table  1.  We  assume  in  this  section,  for 
convenience,  that  A/  .V. 


Step  1 

A'p 

Ng  is  the  number  of  Gauss- 
Legendre  nodes  created. 

Step  2 

0(i/^  A'^  log  .V  +  i/M  log  .V ) 

One  two-dimensional  and  one 
one-dimensional  FFT. 

Step  3 

0(M^) 

Division  of  each  Fourier  coeffi¬ 
cient  FG(m,n)  by  —2nim. 

Table  1:  Time  complexity  of  the  algorithm. 


For  a  given  accuracy  c  the  number  of  Gauss- Legendre  nodes  on  a  single  path  *)/.  accord¬ 
ing  to  Corollary  3.1.  is 

0  ^max  jML.log  .  (44) 

where  M  is  the  highest  frequency  along  7/  and  L  is  the  length  of  -yj  as  described  in  Re¬ 
mark  3.2.  Since  M  =  ,V.  Af  =  .V.  Therefore,  the  number  of  nodes  needed  for  Tj  is: 

O  ^maxj.V  -lirjIl.Lj  ' 7  |)  • 

where  ||  Fj  ||  is  the  length  of  ,  and  Lj  is  the  number  of  the  paths  7^  that  form  Fj .  It  follows 
that  the  total  number  of  nodes.  Ng.  satisfies  the  relation: 

maxj.V  |j 

In  the  worst  case,  i.e.,  when  the  perimeter  is  largest,  there  are  0{N^)  domains  Dj  with 
diameter  0(  1/A’ )  each.  Mg,  then,  becomes 

Ng  =  O  ^max  |  A’  •  A'^  •  A’^  •  log  - .  (47) 
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Choosing  parameters  p  and  according  to  (40),  i.e.,  with 

2i/  =  p  ~  log  (48) 

we  arrive  at  a  bound  for  time  complexity  of  the  algorithm; 

0(p^Ng  +  u'^N^ogN  +  N^)  =  0  ((log^;)  +  (log^  7)  N^logN^  .  (49) 

5  Numerical  Experiments 

We  implemented  the  algorithm  of  section  4  both  with  and  without  the  observation  about  the 
Green’s  theorem  in  §2.2.  The  algorithm  was  implemented  in  FORTRAN  77  and  was  run  on 
SPARCstation  2.  All  internal  calculations  were  performed  in  double-precision  arithmetic, 
but  the  parameters  p  and  u  were  relaxed  for  single-precision  experiments.  Accuracy  was 
measured  in  the  cases  when  all  domains  were  rectangles;  in  these  cases  exact  analytic 
solution  is  available  and  it  was  computed  in  double-precision  for  comparison.  To  assure 
that  integration  along  more  general  curves  does  not  introduce  unexpected  errors,  we  cut 
the  rectangles  into  more  complicated  shapes  (thus  preserving  the  final  answer)  and  verified 
that  the  error  did  not  increase.  This  suggests  that  the  error  is  not  higher  for  the  regions 
for  which  closed  form  solutions  are  not  available. 

In  the  tables  below  we  report  the  results  of  numerical  experiments  using  the  following 
three  algorithms.  In  Algorithm  1  we  make  use  of  the  Green's  theorem  as  described  in 
the  preceding  sections  and,  in  addition,  make  use  of  Remark  4.1  to  integrate  along  the 
horizontal  and  vertical  lines.  Algorithm  2  is  also  based  on  the  Green’s  theorem,  but  without 
preferential  treatment  for  horizontal  or  vertical  fines.  In  Algorithm  3  the  area  integrals  are 
treated  directly  without  recourse  to  the  Green's  theorem. 

W’hen  all  the  domains  {Dj}  are  rectangles,  we  report  run  times  for  a  direct  method, 
which  we  also  use  to  estimate  the  errors  of  the  computations.  When,  for  each  j.  Dj  = 
(a^.hj)  X  (Cj,dj),  the  direct  method  consists  of  calculating  the  Fourier  transform.  0  :  — 

.  by  evaluating  the  Fourier  integrals  analytically: 


O(m.n)  =  ^  dydi 

j  •'“J  •'cj 

(^  —  2-rrim.bj  _  2irima,  \  / 

,  -..L  )■( 


^—2ntndj  _  ^“2jrtncj 


— 2xjn 


The  error  £o&  that  we  report  is  the  maximal  absolute  error  among  all  the  frequencies; 


£00=  max  l/(Tn,n)-0(m,n)l. 

—  M  <m<  M 
-lV<n<iV 


(50) 


W’e  do  not  report  results  for  Algorithm  3  with  N  =  256  because  its  memory  requirements 
were  excessive  for  our  implementation. 

In  the  examples  1  and  2  we  also  include  the  performance  of  the  standard  FFT.  These 
results  are  for  illustrative  purposes  only;  they  can  be  improved  with  the  help  of  .some 
standard  techniques  such  as  Richardson  extrapolation. 
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Example  1  Here  we  compute  the  Fourier  Transform  of  a  single  large  rectangle  (approx 
imately  0.6  by  0.66).  Extrapolating  the  results  for  the  FFT.  we  see  that  single  precision 
accuracy.  E^o  =  10“',  would  result  in  time  complexity  of  more  than  10®  operations,  and 
for  double  precision  we  would  need  more  than  10^®. 

Example  2  In  this  and  the  next  example  we  calculate  the  Fourier  transforms  of  a  piece 
of  a  VLSI  mask  (courtesy  of  Eytan  Barouch).  In  simulations,  these  masks  usucilly  consist 
of  a  large  number  of  simple  geometric  shapes  that  model  the  integrated  circuit.  In  our 
case,  there  are  1215  rectangles  and  424  triangles.  In  this  example  we  compute  the  Fourier 
transform  of  the  rectangles  only  in  order  to  estimate  the  error  of  the  computation.  In  the 
next  example,  we  calculate  the  Fourier  transform  of  all  the  domains,  including  the  triangles. 
The  results  of  this  example  can  be  found  in  tables  4  and  5,  and  those  of  example  3  in  table  6. 


Example  3  In  this  example  we  compute  the  Fourier  Transform  of  1215  small  rectangles 
(the  same  ones  as  in  the  previous  example)  and  424  triangles.  Their  total  area  is  0.183  and 
the  total  perimeter  of  their  union  is  69.3.  The  total  perimeter  of  the  domains  in  this  example 
(and,  consequently,  the  run  times  of  algorithms  1  and  2)  are  sbghtly  smaller  than  in  the 
previous  example  because  some  of  the  newly-added  triangles  share  some  of  their  boundaries 
with  the  rectangles.  (VVe  remove  from  the  computation  the  curves  whose  contributions  can 
be  ezisily  seen  to  cancel.) 

The  following  observations  can  be  made  based  on  the  numerical  experiments  above: 

1.  The  errors  for  all  three  algorithms  do  not  grow  with  N — the  size  of  the  problem. 

2.  All  three  algorithms  have  approximately  the  same  accuracy. 

3.  The  run  times  grow  approximately  like  in  example  1  for  all  three  algorithms,  and 
in  example  2  for  algorithms  I  and  3  (in  agreement  with  (49)),  while  the  run  times  for 
Algorithm  2  in  example  2  appear  to  grow  linearly.  This  apparent  linear  growth  is  the 
result  of  the  large  number  Ng  of  Gauss- Legendre  nodes  generated,  causing  Lagrange 
interpolation  to  dominate  the  computation.  Since  Ng  grows  linearly  with  A',  the  total 
run  times  appear  to  grow  linearly  a.s  well.  The  timings  for  Algorithm  1  in  example  2 
grow  quadratically  because  FFT  dominates  the  computation.  On  the  other  hand,  in 
Algorithm  3  in  example  2  Lagrange  interpolation  dominates  the  computation,  but  A’, 
in  this  case  grows  quadratically.  so  the  overall  growth  remains  quadratic.  The  V'LSI 
mask  of  example  3  cannot  be  bandied  with  the  method  of  Remark  4.1  alone:  some  line 
integrals  need  to  be  computed  numerically.  As  a  result  the  run  times  of  Algorithm  1 
in  example  3  grow  at  an  intermediate  rate. 

4.  Algorithm  2  is  significantly  faster  than  Algorithm  3,  even  when  the  ratio  of  perimeter 
to  area  is  large.  But  Algorithm  1  is  only  twice  as  fast  as  Algorithm  2  on  realistic 
problems  (examples  2  and  3).  On  the  same  realistic  problems  all  three  algorithms  are 
faster  (for  A'  >  64)  than  the  closed  form  solution  even  when  it  is  available.  When  the 
closed  form  solution  is  not  available,  straightforward  calculations  are  prohibitively 
slow:  single  precision  would  require  at  least  150  years  and  double  precision  over 
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Double  Precision:  p 

=  16  and 

u  =  s 

Run  times  (seconds) 

The  error  £oc 

N 

Alg.  1 

Alg.  2 

Alg.  3 

FFT 

Direct 

Alg.  1 

Alg.  2 

Alg.  3 

FFT 

16 

2.3 

129 

0.1 

0.2 

4.8e-15 

6.3e-15 

6.3e-15 

32 

11.1 

671 

0.2 

0.8 

4.6e-15 

4.6e-15 

3.3e-15 

2.3e-02 

64 

52.0 

2755 

0.5 

3.0 

2.0e-15 

2.0e-15 

1.6e-15 

9.2e-03 

128 

208.0 

222.0 

11112 

2.0 

12.2 

l.Oe-15 

l.le-15 

256 

1000.0 

1028.0 

8.7 

48.6 

l.Oe-15 

1.2e-15 

2.5e-03 

Table  2:  Run  times  and  errors  for  example  1;  double  precision. 


Single  Precision;  p 

=  10  and 

u  =  b 

Run  times  (seconds) 

The  error  Eoc 

mm 

Alg.  1 

Alg.  2 

Alg.  3 

FFT 

Direct 

Alg.  1 

Alg.  2 

msmM 

FFT 

16 

0.9 

1.3 

24 

0.1 

0.2 

1.7e-08 

1.5e-08 

4.0e-02 

32 

3.9 

4.7 

95 

0.2 

0.7 

8..5e-09 

8.3e-09 

7.7e-09 

2.3e-02 

64 

18.4 

381 

0.5 

3.1 

5.2e-09 

4.7e-09 

9.2e-03 

128 

83.3 

86.4 

1.538 

2.0 

12.2 

2.0e-09 

5.7e-09 

256 

339.0 

8.7 

48.6 

1..5e-09 

4.4e-09 

Table  3;  Run  times  and  errors  for  example  1;  single  precision. 


Double  Precision:  p  =  16  and  =  8 


Run  times  (seconds) 


The  error  £ 


N 

Alg.  1 

Alg.  2 

Alg.  3 

Direct 

Alg.  1 

Alg.  2 

Alg.  3 

FFT 

16 

6.4 

77 

283 

3.1 

165 

l.le-14 

l.Oe-14 

5.9e-15 

32 

15.8 

156 

887 

3.1 

672 

6.2e-15 

9.4e- 15 

7.7e-15 

64 

55.0 

293 

2277 

3.5 

2707 

5.7e-15 

l.le-14 

5.1e-15 

3.  Ip-02 

128 

213.0 

643 

6562 

4.8 

10874 

7.8e-15 

3.5e-15 

256 

1004.0 

1834 

11.7 

43740 

Table  4:  Run  times  and  errors  for  example  2;  double  precision. 


Single  Precision:  p 

=  10  and 

i/  =  5 

Run  times  (seconds) 

The  error  E,^ 

■■ 

Alg.  1 

Alg.  2 

Alg.  3 

FFT 

Direct 

Alg.  1 

Alg.  2 

16 

2.9 

31 

94 

3.1 

165 

2.2e-08 

3.8e-08 

1.3e-08 

32 

6.0 

48 

190 

3.1 

673 

2.2e-08 

1.8P-08 

64 

19.6 

83 

398 

3.5 

2717 

1.3e-08 

1..3p08 

128 

85.2 

107.5 

4.8 

10920 

9.2e-09 

1.6P-08 

9.0e-09 

256 

.338.0 

562 

BB 

43630 

5..3e-09 

2.7e-08 

1.7e02 

Table  .5:  Run  times  and  errors  for  example  2:  single  precision. 


Double  precision: 
p  =  16  and  u  =  8 

Single  precision: 
p  =  10  and  =  5 

Alg.  1 

Alg.  2 

Alg.  .3 

Alg.  1 

Alg.  2 

Alg.  3 

16 

11 

7.3 

333 

4.2 

29 

32 

23 

148 

989 

8.1 

45 

216 

64 

65 

282 

2451 

22.3 

80 

441 

128 

231 

628 

6972 

89.7 

1 1 5.3 

256 

10.36 

1804 

.347.0 

I 


Table  6:  Run  times  (in  seconds)  for  example  U. 


10'*’  years  regardless  of  the  values  of  .V.  Adjustments  such  as  repeated  Richardson 
extrapolation  will  undoubtedly  improve  these  times,  but  are  unlikely  to  make  them 
competitive. 

6  Conclusions  and  Generalizations 

6.1  Generalizations 

The  entire  algorithm  of  Section  4,  including  Remark  4.1.  readily  generalizes  to  higher  di¬ 
mensions:  we  will  not  describe  the  details  of  the  higher-dintensional  implementation  here. 

We  will  describe  an  extension  of  the  algorithm  to  piecewise  smooth  functions.  If  the 
function  f{x.y)  'm  ( 13)  is  piecewise  smooth,  i.e..  for  each  j,  the  coefficient  Kj  =  A'j(z-.j/)is 
a  smooth  function  of  (z.y),  we  analyze  it  into  a  Fourier  series  with  respect  to  t: 

=  (.'ll) 

m 

where 

l<j{fn){y)^=  /' 

Jo 

We  observe  next  that  all  the  computations  of  1(3. '2  continue  to  hold  if  h'j  is  a  function  of  y 
alone.  That  is.  the  algorithm  of  Section  4  can  be  appbed  to  A'jlrr>)ly)  separately  for  each 
m  and  the  results  added  up  afterwards.  The  time  complexity  of  this  new  algorithm  would 
be  that  of  the  algorithm  of  .Section  4  times  the  width  of  the  widest  Fourier  transform  (in 
the  T-direction  only)  of  the  functions  h'j(x.y). 

6.2  Conclusions 

We  have  presented  an  algorithm  for  evaluation  of  the  Fourier  transform  of  piecewise  constant 
functions  of  two  variables.  The  algorithm  overcomes  the  accuracy  problems  associated  with 
computing  the  Fourier  transform  of  discontinuous  functions:  in  fact,  its  time  complexity 
is  0(log^(  1/c  ).V^  log  .V ),  where  <  is  the  accuracy  and  A'  is  the  size  of  the  problem.  The 
algorithm  is  based  on  the  Lagrange  interpolation  formula  and  the  Creen's  theorem,  which 
are  used  to  preprocess  the  data  before  applying  the  FFT.  It  admits  natural  generalizations 
to  higher  dimensions  and  to  piecewi.se  smooth  functions. 


References 

[1]  .M.  Abramowitz  and  I.  .Stegun,  Handbook  of  Mathematical  Functions,  (Dover,  .New 
York.  1970). 

[2]  E.  Oran  Brigham.  The  Fast  Fourier  Transform  and  its  Applications.  (Prentice  Hall, 
Englewood  Cliffs,  NJ.  19KH). 

[3]  .1.  W.  Cooley  and  J.  W.  Tukey,  An  algorithm  for  the  machine  computation  of  complex 
Fourier  series.  Math.  Comp.  19.  (196-5). 


At, 


* 


♦ 


» 


B 


B 


17 


[4]  G.  Dahlquist  and  A.  Bjorck,  Numerical  Methods,  (Prentire-Hall.  Englewood  Cliffs.  .NJ. 
1974). 

[5]  D.  Gottlieb,  M.  Hussaini,  and  S.  Orszag,  in  Spectral  Methods  for  Partial  Differential 
Equations,  edited  by  R.  G.  Voigt,  D.  Gottlieb,  and  M.  Y.  Hussaini.  (SIAM,  Philadelphia 
PA.  1984). 

[6]  1.  S.  Gradshteyn  and  I.  M,  Ryzhik,  Table  of  Integrals,  Series,  and  Products,  (Academic 
Press,  .New  York,  1980). 

[7]  J.  Stoer  and  R.  Buiirsch,  Introduction  to  Numerical  Analysis,  (Springer  Verlag,  New 
York,  1980). 

[8]  H.  .Joseph  Weaver,  Theory  of  Discrete  and  Continuous  Fourier  Analysis,  (W'iley,  1989). 


